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INTRODUCTION 


The program presented here uses a Monte Carlo technique to 
compute the geometrical factor of a detector and the spatial and 
angular distributions of simulated particles at designated layers 
within the detector. As the program is written, each layer is 
assumed to be rectangular, parallel to every other layer with its 
edges parallel to the edges of every other layer, and centered about 
a catmon axis. See Figure 1. In a previous publication,^ the physical 
significance of the geometrical factor, and of the effective geometrical 
factor as a function of path length, is described. 

The distributions are computed by accumulating data frcm a large 
number of individual, simulated events. Each simulated particle is 
assigned independent, random spatial coordinates (XO, YO) in 
the plane of the uppermost layer of the detector and independent, random 
trajectory angles (THETA, PHI) . A pictorial representation of a 
simulated event in a hypothetical detector configuration is shewn in 
Figure 2. 

This program is written in Fortran IV for the IBM 360 computers and 
is currently being run under the Fortran IV (H) compiler. Seme features 
such as the subroutine ERRSET, RANDU, and REMTIM and the type declara- 
tions REAL *8 and LOGICAL *1 may need to be changed for use with a 
different computer system. 

The operation of the program is described in detail in the next 
section. In the third seertian the use of input data cards and a sample 
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data set are presented. The symbol definition list is given in 
Appendix A. Card images of the Fortran object deck and the graph- 
ing subroutine are listed in Appendix B. A flow diagram of the 
program is presented in Appendix C. 

DESCRIPTION OF THE PROGRAM 

The program begins by finding the CPU and I/O times used in 
the computations for the previous detector configuration by a call 
to the subroutine KEMTIM. Then the initial values for a new data set 
are read. A data set consists of the numerical specifications of a 
particular detector configuration for NPAR number of simulated particles. 
Next, the variables NPLANE and TIME are tested. The program is 
terminated if NPLANE is set equal to 99 on the input data card. If 
insufficient CPU and I/O time remain, TIME = 1 and the program is 
terminated. If execution continues, the remainder of the input data 
narrte for the present configuration are read and the internal variables 
are initialized. The longest possible path length within the detector, 
LPPTH, is determined for later reference. 

In the main loop of the program, the trajectory of a simulated 
particle is generated and tested at successively deeper layers until 
the trajectory passes outside the detector. At the start of the loop, 
after every 3000th simulated particle, the remaining CPU and I/O times 
are again determined. If sufficient time does not remain for both an 
additional 3000 events and the output sunmary for the detector config- 
uration , the program branches to the final computations , resetting 
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NPAR to the ninnber of events presently accumulated. If sufficient 
tiire do es remain, TOPT and TOPW are set equal to the dimensions of 
the uppermost or OTH layer of the detector, TO and WO, The coordi- 
nates of a new simulated particle in the uppermost layer of the 

detector and the zenith and azimuthal trajectory angles are chosen 

2 . ... 

by four separate calls to subroutine RANDU. An isotropic distri- 
bution of incident flux is simulated by the algorithms in the 
program as it is presented in the next section. The change, indicated 
in the caiment cards, will generate a CDS 2 (THETA) distribution of 
incident flux. Advantage is taken of the four- fold symmetry in 
azimuth angle, by varying PHI over the range 0 to tt/2 rather than a 
full 2ir radians. 

After arrays for the graphing output are incremented, the next 
major, nested loop of the program is entered. The detector configur- 
ation is treated as a series of boxes with rectangular tops and 
bottcrrs and trapizoidal sides. These boxes, or levels, are numbered 
by the DO variable K the same as the number of the bottom layer or 
plane. The surface through which the particle exits the detector 
level is determined. If the particle has exited the detector, data 
for that particle are stored, and the program returns to generate a ' 
new simulated particle. If the particle has passed through the bottom 
of the present detector level into the next level, the parameters are 
set to the values for the next detector level and the loop is repeated. 

In the last detector level (K = NPLANE) the path length of the trajectory 
within that level is ccnputed. The distribution of path lengths and 
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the effective geometrical factor as a function of path length are 
output at the carpletion of the data set. 

After NPAR simulated particles have been generated and their 
trajectories through the detector levels determined, the output 
sunmary for the detector configuration is prepared. Path length 
data are computed and Written. If the detector has the form of a 
rectangular parallepiped, the geometrical factor is computed 
analytically and written. Finally the distributions of simulated 
particles by path lengths, by spatial coordinates at specified detector 
layers, and by zenith angle at specified detector layers are printed 
numerically and as printer plots through the use of the graphing sub- 
routine. Following the output execution, the program returns to the 
beginning. 

DESCRIPTION OF THE INPUT DATA 

One data set, including a heading, the detector dimensions, and 
other numerical parameters, is specified for each detector configur- 
ation. All dimensions are in centimeters. The data cards to be read 
for each data set are as follows. 

Card 1 a heading. Any combination of 80 characters and blanks an the 
first card will be read and written out at the start of the output. 
Card 2 detector and particle data. The input data are read in the 
following order: the number of detector layers after the first, 

NPLANE, in format 12; the X dimension of the 0th plane, TO, in 
format F10.4; the Y dimension of the 0th plane, WO, in format F10.4; 
the number of simulated particles to be generated, NPAR, in format 
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110; the initial value for the randan integer, IY, in format 15; 
and the numbers of those detector layers for which distribution 
data are required, up to 10 plane numbers in ascending order, in 
format 1012. 

Card 3 and subsequent cards up to a total of NPLANE +2. One card is 

prepared for each detector layer after the topmost, giving in 

order the X dimension of the layer, T(*) ; the Y dimension of the 

layer, W(*) ; and the distance along the Z axis betoeen the present 

and the previous layers, H(*) , in format 3F10.4. 

Note that after the data set for the last detector configuration, two 

more data cards are needed. The first card, a fake heading, is followed 

by another card with the digits 99 in the first too columns. This sets 

NPLANE to 99 which terminates the program. A sample data set follows: 

THIS DETECTOR IS A THREE -EIEMENT COSMIC RAY TELESCOPE 
2 94.6 94.6 100000 37 1 2 

76.5 76.5 18.2 

76.5 76.5 57.3 

THIS IS A FAKE HEADING CARD 
99 

This data set specifies a detector configuration with three layers. 

The top layer is 94.6 x 94.6cm 2 . A total of 100000 simulated particles 
will be generated using the random number generator initialized with 
the integer, 37. Spatial and angular distribution data will be generated 
for the first and second layers as well as for the topmost. The first 
layer following the top has dimensions 76.5 x 76.5 on^ and is spaced 18.2 cm 
below the top. The second layer has dimensions 76.5 x 76.5 cm and is 
spaced 57.3 an below the first. The program is terminated after finishing 
the output surrmary for this detector configuration. 
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SYMBOL DEFINITION LIST 



AND ( 30 ) I int| the number of particles that pass through the last plane of 
the detector divided into 30 three-degree bins in THETA. 

ANGD 15(4/4,30/11) |int| the distribution of particles by position and by 
angle with respect to the Z axis at 11 different planes. _ The 
subscripts represent: the displacement from the X axis divided into 
4 equal bins; the displacement from the Y axis divided into 4 equal 
bins; the angle made with the Z axis, THFTA, divided into. 30 three- 
Hegree bins; the Nth plane number given by GRAPM(N), N=1 is always 
the 0th plane. 

AM0DK4) |int| every other 9-degree bin, 3 three-degree bins in ANGDIS, 

starting at 9 degrees. After each set of its values are written out/ 
it is filled with the next bin corresponding to the 2nd subscript of 
ANGDIS/ rotating 1/ 2, 3, 4. 

ANGD2 (4) lintl the same as ANGD1 except it represents the alternate 
9-degree bin starting at 18 degrees. 

ANGLE Ireall set to every integer multiple of 3., up to 90., for reference 
in output. 

ANGT lintl used as the 3rd subscript to ANGDIS to divide the angle made 
with the Z axis, THETA, into 30 three-degree bins. 

ANGX lintl used as the 1st subscript to ANGDIS to divide a particle’s 
displacement from the X axis, plus or minus, into 4 equal bins. 

ANGY lintl used as the 2nd subscript to ANGDIS to divide a particle's 
displacement from the Y axis, plus or minus, into 4 equal bins. 

AVEPTH Ireall the average path length of all particles that enter the last 
detector level . 

CGF Ireall the geometrical factor which can be computed analytically for 
the special case that the detector shape is a rectangular 
paral lei epl ped. 

CPU I Ireall cos (PHI) 

CPSD Ireall the particle counts per square degree. 

CTHETA Ireall cos(THFTA) 

DEG1 Ireall set to everv other multiple of 9., starting at 9., up to 81., 
for output reference. 

DEG2 Ireall set to everv other multiple of 9., starting at 18., up to 90., 
for output reference. 
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EFFGF Ireall the effective geometrical factor for oath lengths greater 
than J/100.*LPPTH, 0=1,100. 

ERRSET an I P>M system subroutine that suppresses the output which would 
otherwise be printed for error number IHC210. It is possible to 
receive this error when calculating PTHLNG in PATH1 or PATH2. It 
has no significnace in this program because the standard fixup sets 
the result to zero, which is ignored. See * ' I BM System/360 Operating 
System, FORTRAN IV (G & H) Programmers' Guide'', (G028-6217-2) , pp. 
124-125. 

EXITPT(99,3) I inti gives the plane through which the particle leaves the 

detector. The subscripts indicate the detector level and which surface 
within the level is the exit plane of the particle. . The subscripts 
represent; the detector level where the particle exits; the plane 
within the level that the particle exits from as follows: 

(1) if this is the last detector level and the particle passes 
through the bottom detector layer. If this is any other level, 
the particle does not exit the detector by passing through the 
lower layer. 

(2) the particle pass out of either side between the detector layers. 

(3) the particle passes through the front or hack between the 
detector layers. 

G a subroutine that, will do printer plots of up to 100 points. The 
arguments that must he passed to it are as follows; 

TYPE | i n 1 1 with a value greater than 0 as follows, 

(1) a 1 i near plot . 

(2) a semi-log plot In X, and a linear plot in Y. 

(3) a semi-log plot, in Y, and a linear plot in X. 

(4) a log-log plot. 

(5 or greater) the subroutine will decide what kind of plot to 
make. 5 Is used in all calls In this program. 

HEAD(5) |real*8| should contain the heading for the plot. 

A(N) |real |, N between 1 and 100, contains the X coordinates of the 
data. 

R(N) ! rea 1 | , N between 1 and 100, contains the Y coordinates of the 
data . 

SIZE ( i n t ! contains the dimensions of A and R. 

GAND ( 30 ) |real| a variable that transmits the values of AND to subroutine G. 

GANG(30) |real| a variable that transmits the values of ANGLE to subroutine 

G. 

GCPS ( 30 ) |real| a variable that transmits the values of CPSD to subroutine 
G • 

GEFF(100) |real| a variable that transmits the values of EFFGF to 
subroutine G. 
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GF |real | the geometrical factor of this detector for particles passing 
through every detector layer. 

GPBI(IOO) |real| a variable that transmits the values of PBIN to subroutine 
G. 

GPOL ( 30 ) | real | a variable that transmits the values of POLANG to 
subroutine G. 

GRAPH ( 11 ) | int| contains the plane numbers for the last subscripts of 

ANGD IS, XDIS, and YD I S . GRAPH ( 1) always has the value of 0, GRAPH 
must be filled left, to right In ascending order. 

GTAN ( 30 ) | real | a variable that transmits the value of TANG to subroutine 
G. 

GTOK30) |real| a variable that transmits the value of TOTC to subroutine 

G • 

GXB ( 50 ) I real I a variable that transmits the values of XC to subroutine G. 

GXD 1(50) Ireall a variable that transmits the values of XDIS to subroutine 
G. 

GYB( 50) Ireall a variable that transmits the values of YB to subroutine G. 

GYD 1(50) Ireall a variable that transmits the values of YDIS to subroutine 

G. 

H(99 ) Ireall the Z dimension of the Nth detector level. 

HEAD ( 20 ) Ireall reads and writes the data set heading. 

HO Ireall is the total height of the detector in the Z dimension in the 
case that the detector has shape of a rectangular parallelepiped. 

I I i nt I a DO var i able. 

I CPU | int| the first argument to REMTIM, it returns the seconds of CPU 
time remaining. 

ICPUS I int| saves the previous value of ICPU. 

I 10 | inti the 2nd argument to REMTIM, it returns the seconds of I/O time 
remai ning. 

NOS lint) saves the previous value of NO. 

IX | inti is the 1st argument to REMTIM, it should contain the previous IY. 
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IY | Inti the 2nd argument to RANDU, it Is a random integer that is returned 
by RANDU. 

11 | i nt I a DO var i able. 

12 1 1 nt: I initial value for II in a DO loop. 

13 I inti the maximum value for II in a DO loop. 

J I int | a DO variable. 

JS | I nt | saves the current value of J for later use. 

K | i nt | a DO var i abl e. 

K5 | I nt I saves the current value of K for later use. 

L | i nt | a DO var i able. 

LAPTH |real| the longest actual path length of any particle that enters 
the last detector level. 

LP |real| used to calculate LAPTH. 

LPPTH |real| the longest possible path length in the last detector level. 

L5 | i nt | saves the current value of L for later use. 

M | int | a DO variable. 

M | int| a DO variable. 

fJPAR | Int| the number of particles to he simulated in this data set. 

NPLANF | int I the number of planes after the Oth. NPI..ANF = 99 stops the 
program. 

PATH1 Ireall a statement function subprogram which is used to calculate 
PTHLNO . 

PATH2 |real| a statement function subprogram which is usnd to calculate 
PTHLNO . 

PBIM |real| the longest possible path length, I.PPTH, divided into 100 
eaua 1 bins. 

PHI lreal| the azimuthal angle. 

POLAR I i nt | set to even multiples of 3, from 3, up to 90, for output 
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reference . 


POS | inti used to determine a particles position with respect to the edees 
of the bottom of the present detector level: 

(1) Inside the edges of the plane. 

(2) to the right or left of the plane, 

(3) to the front or back or the plane, 

(4) to both the front or back and right or left of the plane, that 
is, off a corner. 

PTHI.NG | real | the length of the path of a particle In the last detector 
level . 

GC |real| an Internal variable for calculating X and Y. 

Ql I inti the subscript to AND. It divides THETA into 30 three-degree bins. 

QPTH |real| an internal variable used to calculate PTHI.NG. 

QPTH1 (real | an internal variable used to calculate PTHING. 

QPTH11 (real | an internal variable used to calculate PTHING. 

QPTH2 [real | an Internal variable used to calculate PTHLNG. 

QPTH22 | real | an Internal variable used to calculate PTHLNG. 

QS |real| an Internal variable used to accumulate the values of AND In 
order to suppress zeros. 

QX |tnt| used as the 1st subscript to XDIS, It divides TO Into 50 eaual 
bins . 

QY 1 1 nt I used- as the 1st subscript to YDIS, It divides WO into 50 equal 
bins. 

RANDU a subroutine which Is used to generate random numbers. The arguments 
are as follows: 

IX I Inti should be set to the previous IY, RANDU uses IX to generate 
another IY and YFL. 

IY | i nt | is the random Integer returned. It Is read Into the 

program as Input data. Its Initial value Is any odd Integer of 
9 or less digits. 

YFL Ireall is a random number between 0.0 and 1.0 which is used to 
generate particle parameters. 

REM I Inti used to count 3000 particles between calls to REMTIM. 

SIZE I Inti the last argument to G which contains the dimension of GAND 
and GTAN . 
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5PH I | real I sJ n(PIII ) 

SQX Ireall an internal variable used to calculate CGF. 

5QY |real| an internal variable used to calculate CGF. 

STHETA Ireall sin (THETA) 

SUBHK5) |real*8l contains the heading for the 1st call to G, 

"EFFGF VS. PTHLNG". 

5UBH2(5) | real * 8 1 contains the heading for the 2nd call to G, 

"TOTAL COUNTS VS. PATH-LENGTH". 

5UBH3( 5) | real *8 1 contains the heading for the 3rd call to G, 

"COUNTS OUT BOTTOM VS. POLAR ANGLE". 

SUBH4( 5) I real *8 1 contains the heading for the 4th call to G, 

' 'CPSD VS. POLAR ANGLE" . 

SUBH5 ( 5 ) I real *8 1 contains the heading for the 5th call to G, 

"NO. OF PARTICLES VS. XBIN". 

SUBH6(5) | real * 8 1 contains the heading for the 6th call to G, 

"NO. OF PARTICLES VS. YBIN". 

SUBH7(5) |real*8| contains the heading for the 7th call to G, 

"TOTAL PARTICLES VS. ANGLE". 

SUMC(IOO) I inti the number of particles which have a path length 
in the last detector of N/100.*LPPTH or greater. 

T(99) Ireall the X dimension of the Nth plane. 

TANG 1 1 nt I the total number of particles that are within a specific 
3-degree angular bin at a given plane. 

TGR(ll) |real| the X dimension of the Nth plane in GRAPH(N). 

THETA |real| the zenith angle of a simulated particle trajectory. 

TIME I i nt | a flag/ when set to one, the program stops. 

TOPT | real | the X dimension of the top of the detector level that the 
particle is presently passing through. 

TOPW |real| the Y dimension of the top of the detector level that the 
particle is presently passing through. 
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TOPX |real| the X coordinate of a particle at the top to the present 
detector level . 

TOPY |real| the Y coordinate of a particle at the top of the present 
detector 1 evel . 

1010(4,100) | l nt | used to determine the path length and exit place of 
the particles from the detector. The subscripts represent: the 
exit pi ane, 

(1) bottom of the detector, 

(2) side of the detector, 

(3) front or back of the detector, 

(4) the total of 1, 2, and 3; 

path length divided into 100 equal bins. 

TOTPTH |real| the sum of all of the path lengths of all the particles 
that enter the last detector level. 

TTANG |lnt| the total of all of the particles In a given plane. TTANG 
is used to test against ZF.RO In order to suppress zeros. 

TO |real| the X dimension of the 0th plane. 

W(99) lreal| the Y dimension of the Nth plane. 

WGR(ll) |real[ the Y dimension of the Nth plane in GRAPM(N). 

W I CPUS |real| used to write out, in minutes, the CPU time required by 
the previous data set. 

WII05 (real | used to write out, in minutes, the 1/0 time required by 
the previous data set. 

WLPPTH |real| used to write out the longest possible path length. 

W0 |real| the Y dimension of the 0th plane. 

X(99) |real| the coordinate of a particle in the X dimension at the Nth 
pi ane. 

XB |reall set to multiples of the X dimension of a plane divided by 
100, for output reference. 

XDIS(50,11) |lnt| the distribution of particles along the X axis, that 
is in strips parellel to the Y axis, the subscripts represent: 
the displacement from the X axis, plus or minus, divided into 50 
equal bins; the Nth plane in GRAPH(N). 

X0 | real | the X coordinate of a particle at the 0th plane. 
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Y ( 99 ) Ireall the coordinate of a particle In the Y dimension at the Nth 
plane. 

YD (real I set to a multiple of the X dimension of a plane divided by 
100, for output reference. 

YDIS(50,11) lintl the distribution of particles along: the Y axis, that is 
in strips parellel to the X axis, the subscripts represent: the 
displacement from the Y axis, plus or minus, divided into 50 equal 
bins; the Nth plane in GRAPH(N). 

YFL [real | an argument to RANDU, it Is returned with a random real 
number between 0.0 and 1.0. 

Y0 |real| the Y coordinate of a random particle at the 0th plane. 

ZERO lintl used to find the number of particles at each plane in GRAPH(N) 
and to suppress zeros in the output. 
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LISTING OF THE 


GEOMETRICAL FACTOR PROGRAM 
AND 


GRAPHING SUBROUTINE 



RFAL T(99),U(99),H(9q),X(99),Y(99),HFAn(20),TRR(ll),WGR(ll), 

AGPRI (mn),GFFF(100),RPnL(30),GCPS(30),GANn(30),GTOT(100),RTAN(30), 
P.GANR(30),RXR(50),RYR(50),RXni (50),RYni ( 50) , LAPTH, LPPTH, LP 
INTEGER EXITPT(99.3),POLAR,AND(30),SUMC(100),TOTC(4, 100),OX,OY, 
AXDIS<50,11 ), YD I S( 50, 11 ) , GRAPH ( 11 ) , ANGX, AMRY, ANGT, ANGD I 5(4,4,30,11) 
R, ZERO, TANG, ANRD1 ( 4 ) , ANGD2 ( 4 ) , TTANG, T I ME, POS, 0 I , RFM, S I ZE 
COMMENT (SUBH1-7) CONTAIN THE HEADINGS FOR TRANSMISSION TO SUBROUTINE 
C (G) THEY ARE REAL* 8 FOR CONVINENCE ONLY. 

REAl.*8 SUBH1 ( 5 ) , SUBH2 ( 5 ) , SIJBH3 ( 5 ) , SUBH4 ( 5 ) , SUBH5 < 5 ) , SUBH6 ( 5 > , 
ASI.IRN7C5) 

DATA SIIBH1, SIJRH2, SUBH3, SMRH4, SIIRH5, SUBH6, SUBH7/'EFFRF VS','. PTHLN 
AG', 3*' ', 'TOTAL CO','IINTS VS.',' PATH-LE ' , ' NOTH 

R ' , 'COUNTS 0','lJT BOTTO ' , 'M VS. PO','I.AR ANRL ' , ' F ','CPSD 

C VS.',' POLAR A','NGLE ',2*' ','NO. OF P ' , 'ARTICLES ' , ' V 

DS. X B * , * I N ',' ','NO. OF P ' , ' ART I CLES ' . ' VS. Y B','IN 

E ', 'TOTAL PA ' , ' RT I CLES '.'VS. ANRL ','F '.' 

F '/ 

COMMENT (PATH1, PATH2) ARF STATMENT FUNCTION SUBPROGRAMS WHICH 
C CALCULATE A COEFFICIENT FOR THE PATH LENGTH OF THE PARTICLES. 

PATH1 (A, B,C.D)=(A-2.*R)/(A-C+2.*(D-B)) 

PATH2 ( A, B,C,D)=(2. *A+B )/(B-C+2.*(A-D)) 

COMMENT (T I ME ) IS A FLAG, WHEN SET TO ONE THE PROGRAM STOPS 
TIMF=0 

COMMENT THIS DO Al I OWS THE PROGRAM TO OPPFRATE ON AS MANY AS 100000 
C DATA SETS. 

DO 10 1=1,100000 

COMMENT (RFMTIM) IS AN IBM SUBROUTINF THAT RFTURNS THE SECONDS 
C REMAINING OF CPU AND I/O TIME AS INTEGERS. 

CALL REMT IM( I CPU. I 10) 

I F ( I .FO.DRO TO i.01 

COMMENT THIS DETERMINES THF TIMF REOll I RED FOR THF PRFVIOIIS DATA SET 

C AND WRITFS THIS OUT AS MINUTES OF CPU AND I/O TIME. 

W I C PUS - ( I C PUS - I C PIJ ) / 6 0 . 

Wl IOS=(l IOS-1 I0)/60. 

WRITE (6, 1000 )WI CPUS, Wl I OS 
101 ICPUS-ICPU 
I IOS= I 10 

COMMENT SET (ANGDIS, GRAPH) TO ZERO. 

DO 11 J-1,11 
DO 12 K=l, 30 
DO 12 L=l, 4 
DO 12 M= 1 , 4 
12 ANGD IS(M,L,K,J)=0 
11 GRAPH (J ) = 0 

COMMENT READ THE DATA SET HEADING, THE NUMBER OF PLANES AFTER THE OTH, 

C LENGTH OF THE OTH PLANE IN THE X DIMENSION, THE LENGTH OF THE 

C OTH PLANE IN THE Y DIMENSION, THE NUMBER OF PARTICLES, THE 

C INITIAL (I Y), THE PLANES FOR WHICH DISTRUBUTION DATA IS 

C REQUESTED IN ASSENDING ORDER. THE FORMAT IS, (20A4/ 1 2, 2F10. 4, 

C 110,15,1012). 
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READ (5*1001) HEAD, NPLANE , TO, WO, N PAR, I Y , (GRAPH (J ), J = 2, 11 ) 

COMMENT TEST TO SEE IF THERE ARE NO MORE DATA SETS, IF (NPLANE )=99, 

C OR IF (TIME) HAS BEEN SET TO ONE. IF SO BRANCH TO THE STOP 

C STATEMENT. 

I F ( (N PLANE . EO,. 99 ) . OR . (T IME . EO. 1 ) )GO TO 102 
WRITE (6, 1002)1, HEAD, NPI. AN E,N PAR, IY,T0,VJ0 
COMMENT INITI.IZE (OS, LAPTH, HO . XD I S . YD I S, TOPT, EX I TPT, GAND, GANG , GCPS 
C GPOL,GTAN,AND,TOTC,GEFF,GPBI ,GTOT,SUMC) TO ZERO 

QS = 0 . 

LAPTH=0. 

H0 = 0 . 

DO 13 J-1,11 
DO 13 K=l, 50 
XD I S ( K, J ) =0 

13 YD I S ( K, J ) = 0 
TOTPTH=0 . 

DO 14 J=l, 99 
DO 14 K=l,3 

14 EX I TPT ( J , K) = 0 
DO 30 J=l, 30 
GAND ( J ) = 0 . 

GANG ( J ) s 0 . 

GCPS(J ) = 0. 

GPOL ( J ) =0 . 

GTAN(J )=0. 

30 AND ( J )=0 

DO 15 J=l,100 
DO 10 K-1,4 
10 TOTC(K,J)=0 
GEFF ( J ) e 0 . 

GPB I ( J ) =0 . 

GTOT ( J ) =0 . 

15 SUMC ( J ) =0 

COMMENT READ AND WRITE (T,W,H), THE X, Y, AND Z DIMENSIONS OF THE 

C DETECTOR RESPECTIVELY. FOR EACH PLANE AFTER THE OTH, (H) IS 

C THE DISTANCE FROM THE PREVIOUS PLANE. THE FORMAT IS, (F10.4). 

DO 17 0-1, NPLANE 
READ ( 5, 1003)T(J),W(J),H(J) 

17 WRITE(G,1004)J,T(J),W(J),H(J) 

COMMENT SET THE VALUES OF (TGR, WGR) EQUAL TO (T,W) AT EACH PLANE 
C LISTED IN (GRAPH), AS WELL AS TO THE OTH PLANE. 

TGR( 1)=T0 
WGR(1)=W0 
DO 18 0-2,11 

I F ( G RA PH ( 0 ) . E 0 . 0 ) GO TO 100 
TGR(0)=T (GRAPH (0) ) 

18 WGR( .1 ) = W(GRAPH ( J ) ) 

COMMENT DETERMINE 100. /(THE LONGEST POSSIBLE PATH LENGTH), ( LPPTH ) 

100 I F (NPLANE . EO . l)GO TO 104 

I.PPTH-100 . / SORT ( ( (T (NPLANE ) /2+T ( NPLANE-1 ) / 2 ) **2 ) + 


R-2 



A((W(NPI.ANF)/2+W(NPLANF-l)/2)**2)+H(NPLANF)**2) 

LP=100 . /SORT (W(NPLANE-1 ) **2+T ( NPLANE-1 ) **2) 

104 LPPTH*ilOO . / SORT (■( (T ( 1 ) / 2+T0/2 > **2) + ( (W( 1) / 2+W0/2 )**2 )+H ( 1) **2) 

I P =100 . / SORT (T0*T0+W0*W0) 

105 I F ( LPPTH . GT. LP) LPPTH=LP 
REM=0 

COMMENT THIS DO LOOP FINDS A RANDOM PARTICLE AND SENDS IT THRODOH 

C EACH DETECTOR LEVEL WHILE COLLECTING INFORMATION ABOUT IT 

DO 19 J = l, NPAR 

COMMENT AT EVERY 3000TH PARTICLE CALL (RFMTIM) AND TEST TO SEE IF THERE 
C IS ENOUGH TIME FOR ANOTHFR 3000 PARTICLES AND TO FINISH THE 

C OUTPUT SECTION OF THIS PROGRAM. IF NOT GO DIRECTLY TO THE 

C OUTPUT SECTION USING ONLY AS MANY PARTICLES AS HAVE ALREADY 

C BEEN PROCESSED. THE INDICATED TIMES, 14 AND 5, ARE APPROPIATF 

C TO MOST DATA SETS ON AN IBM 360/91 AND SHOULD BF MODIFIED FOR 

C DIFFERENT MACHINES. 

■REM-RFM+1 

IF(REM.NE. 3000)00 TO 103 
CALL REMTIMC ICPU, I 10) 

!F(( ICPU. LE.14) .OR. ( I IO.LF.5))GO TO 106 
COMMFNT SF.T THE VALUES OF (TOPW,TOPT, I X ) ; INITLIZE ( PTHLNG ) 

REM=0 

1.03 TOPW-WO 
TOPT=TO 
PTHLNG=0. 


I X = I Y 

CALL RANDU( IX, I Y, YFL) 

COMMENT FIND A RANDOM (YO) ANYWHERE WITHIN THE OTH PLANE 
Y0=W0* ( YFL- . 5 ) 

I X= I Y 

CALL RANDU( IX, I Y,YFL) 

COMMENT FIND A RANDOM (XO) ANYWHERE WITHIN THE OTH PLANE 
X0=T0*(YFL-. 5) 

I X= I Y 

CALL RANDU( IX, I Y,YFL) 

COMMENT FIND A RANDOM (THETA, STHETA,CTHETA) FOR AN ISOTROPIC 
C D I STURBUT I ON OF INCIDENT FLUX. THE FOLLOWING EXPRESSIONS FOR 

C (STHETA, CTHETA) , WHICH ARE LISTED AS COMMENTS, MAY REPLACE 

C THOSE LISTED AS FORTRAN, TO OBTAIN A COS-SQ DISTRIBUTION OF 

C INCIDENT FLUX. 

C STHETA® SQRT (.999999- . 99999+SQRT (YFL) ) 

C CTHETA=SQRT (SORT (.000001+. 99999 *YFL)) 

CTHETA=SQRT( . 999999-. 99999*YFL) 

STHETA=SQRT( .000001+. 99999+YFL) 

THETA=ARS IN (STHETA) 

I X = I Y 

CALL RANDU( IX, IY,YFL) 

COMMENT FIND A RANDOM (PH I , SPH I , CPH I ) FOR ANGLES UP TO 90. DEGREES 
PH 1 = 1. 570796 *YFL+. 000001 
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SPH I =S I N ( PH I ) 

CPU I =COS (PH I ) 

COMMENT DETERMINE IF THERE ARE ANY VALUES IN (GRAPH) , IF NUT SKIP THIS 

C SECTION WHERE DISTRIBUTION DATA IS COLLECTED FOR THE OTH PLANE 

I F (GRAPH ( 2 ) . EO . 0 )GO TO 107 

COMMENT (ANGT, ANGY, ANGX) ARE USED AS SUBSCRIPTS, (ANGT) DIVIDES THETA 
C INTO 30 EQUAL PARTS, 19. 09859=NUMBER OF DEGREES IN A RADIAN/3. 

C (ANGY,ANGX) EACH PUT ALL OF THE PARTICLES INTO 4 EOUAL BINS ON 

C THE BASIS OF THEIR DISTANCES FROM THEIR RESPECTIVE AXES. 

ANGT=THETA*19 . 09859+1 
ANGY»ABS (Y0)*8/W0+l 
ANGX=ABS(X0)*8/T0+1 

COMMENT (ANGDIS) ACCUMULATES THE ANGULAR DISTURBUTION BY X AND Y BIN 
C AT 11 DIFFERENT PLANES, IT INCLUDES ALL PARTICLES THAT PASS 

C THROUGH EACH PLANE AND, IN THIS CASE, THROUGH THE OTH PLANE. 

ANGD I S (ANGX, ANGY, ANGT , 1)=ANGD I S (ANGX, ANGY, ANGT, 1)+1 
COMMENT (QX,QY) ARE USED AS SUBSCRIPTS TO (XD I S, YD I S ) . THEY DIVIDE 

C THEIR RESPECTIVE AXES INTO 50 BINS NUMBERED FROM THE AXES OUT. 

C SO THAT THE BINS ARE STRIPS, RUNNING PERPENDICULAR TO THEIR 

C RESPECTIVE AXES. 

QX=ABS(X0)*100/T0+1 
QY“ABS ( YO ) * 100/W0+ 1 
XDIS(QX,l)=XDIS(OX,l)+l 
YDIS(OY,l)«:YDIS(QY,l) + l 

COMMENT IN THE LOOP STARTING AT '00 20' FACH LEVEL IS TREATED AS EVERY 

C OTHER, BECAUSE THF OTH PLANE HAS ITS DIMENSIONS GIVEN IN THE 

C VARIABLES (TO, WO), NEW VARIABLES ARE TEMPORARLY USED. THEY 

C ARE (TOPX,TOPY,TOPT,TOPW) WHICH RESPFCTIVFt.Y STAND FOR THE X 

C CORD I NATE OF THE PART I Cl E AT THF TOP OF THE LEVEL, THE Y 

r - CORD I NATE , THE X DIMENSION, AND THE Y DIMENSION, 

107 TOPX = X 0 
TOPY=Y 0 

COMMENT THIS LOOP IS THE HEART OF THE PROGRAM, IT SENDS THE PART I C LF 
C THROUGH FACH DETECTOR I.FVEI, WHII.F ACCUMIILAT I NG DATA. 

DO 20 K=l, NPI.ANE 
COMMENT SAVE THE VALUE QF ( k ) 

KS = K 

OC=STHETA/CTHFTA*H(K) 

COMMENT FIND THF COORDINATES OF THF PART I OLE IN PLANE K. 

X ( K) “TOPX+CPH I *OC 
Y(K) =TOPY+SPH I *QC 

COMMFNT DETERMINE THE POSITION OF THE PARTICLE RELATIVF TO THE EDGES 
C OF LAYER K. 

C (POS)-l, INSIDE THE EDGES OF THE PLANE 

C ( POS ) =2, TO THE RIGHT OR LEFT OF THE PLANE 

C ( POS ) =3, TO THE FRONT OR BACK OF THE PLANE 

C ( POS ) =4, TO ROTH THE RIGHT OR LEFT AND THE FRONT OR BACK OF 

C THF PLANE 

P05 = 4 

I F ( ABS ( X ( K ) ) . LE.T(K)/2.)POS=POS-2 
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IF(AP.S(Y(K)) . LE.W(K)/2. )P0S = P0S-1 
GO 10(108,109,110,111), POS 

COMMENT SEARCH FOR A VALUE IN (GRAPH) THAT MATCHES THIS Pl.ANF , 

C IF THERE IS ONE SAVF THE SUBSCRIPT OF (GRAPH) THAT IT 

C MATCHES. 

108 DO 21 L=2, 11 
LS-L 

! F (GRAPH ( L ). EQ . K)GO TO 112 
21 CONTINUE 
GO TO 113 

COMMFNT ACCUMULATE (XO I S, YD I S, ANGD I S ) AS IN THE PREVIOUS CASE 

112 QX=ABS (X (K))*100/T(K) + 1 
QY=ABS(Y(K))*100/W(K)+1 

xd i s ( ox , ls ) = xn I S ( OX , LS ) + 1 

YD I S (OY, LS ) =YD IS(0Y,LS) + 1 
ANGX=ARS (X ( K) )*8/T(K)+l 
ANGY=ABS(Y(K) )*8/W(K)+l 

ANGD I S (ANGX, ANGY, ANGT, LS ) =ANGD I S ( ANGX, ANGY, ANGT, LS ) + 1 
COMMENT IF THERE ARE MORF PLANES. REINITIALIZE TO START ANOTHER LEVEL. 
C IF NOT, COMPUTE THE PATH LENGTH, ADD THAT TO THE TOTAL PATH, 

C LENGTH, AND FIND THE CORRECT BIN IN WHICH TO INCREMENT (AND). 

113 I F ( K . NE . NPLANE )GO TO 114 
PTHLNG=H ( K)/CTHFTA 
TOTPTH=TOTPTH+PTHI NG 
OI=THETA*19. 09859+1 
AND(O.I )=AND(0 I ) + l 

GO TO 123 

COMMENT IF THIS IS NOT THE LAST DFTFCTOR LEVEL, INCREMENT (EXITPT) 

C WHICH TELLS WHF.RF THE PARTICLES LEFT THE DETECTOR, THFN GET 

C ANOTHER PARTICLE. 

109 I F (K . NF . NPLANE )GO TO 19 

COMMENT FIND THE PATH LENGTH BY TAKING THE LEAST POSITIVE VALUE OF 

C (PTHLNG.QPTH) . 

PTHLNG=PATHl(TOPW,TOPY,W(K) , Y(K) ) 

QPTH a PATH2 (TOPY , TOPW, W( K) , Y( K) ) 

GO TO 124 

110 I F(K.NE. NPLANE )GO TO 19 

COMMENT FIND THE PATH LENGTH BY TAKING THE LEAST POSITIVE VALUE OF 
C (PTHLNG,QPTH) . 

PTH LNG= PATH 1 (TOPT,TOPX,T (K),X(K)) 

QPTH=PATH2 (TOPX ,TOPT, T (K),X(K)) 

124 IF( (PTHLNG.LT . 0. ) . OR. ( (QPTH . LT. PTHLNG) . AND . ( 0. . LE .QPTH) ) ) 
APTHLNG=QPTH 
GO TO 115 

COMMENT FIND THE PATH LENGTH BY TAKING THE LEAST, POSITIVE VALUE 
C OF ( QPTH 1, QPTH 11, QPTH 2, QPTH 22 ) 

111 QPTH 1= PATH 1 (TOPT, TOPX, T (K),X(K)) 

QPTH11=PATH2 (TOPX,TOPT,T (K),X(K)) 

QPTH 2= PATH 1 (TOPW, TOPY, W( K) , Y( K) ) 

QPTH22=PATH2 (TOPY,TOPW, W( K) , Y( K) ) 
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COMMENT IT IS POSSIBLE THAT WHEN CALCULATING IN ( PATH1, PATH2) A 
C DIVISION CY ZERO WILL BE ATTEMPTED. THIS IS POSSIBLE BECAUSE 

C A PARTICLE MAY BF. PARALIEI. TO ONE OF THF SIDFS OF THE DETECTOR 

C TO WITHIN THF ACCURACY OF THF MACHINE. THIS HAS NO DETRIMENTAL 

C FFFECT ON THE RESULTS. AND A CALL TO (ERRSET), CANCELS THE 

C ERROR MESSAGE. EXECPT FOR A STATEMENT AT THE END OF THE OUTPUT 

CALL ERRSETC210, 256,-1,1) 

IF ((OPTHl.LE.O.).OR. ((QPTH11.LT. QPTH1). AND. (0. .LT.QPTH11))) 
AQPTH1=QPTH11 

IF((QPTH2.LE.O. ).OR. ( (QPTH22. LT.QPTH2) . AND. (0. . LT . QPTH22 ) ) ) 
AQPTH2=QPTH22 

COMMENT WHILE FINDING THE LEAST OF (QPTH1, QPTH2 ) , ALSO SET (POS) TO 

C THE APPROP I ATE VALUE TO REFLECT WHERE THE PARTICLE ACTUALLY 

C LEFT THE DETECTOR. 

I F (QPTH1 . LT . QPTH2 )GO TO 116 
I F ( K . EQ. NPLANE ) PTHLNG-QPTH2 
POS = 2 
GO TO 115 

116 I F(K.EQ. NPLANE )PTHLNG-QPTH1 
POS = 3 
GO TO 115 

COMMENT RESET (TOPX,TOPY,TOPW,TOPT) TO VALUES FOR THE NEXT DETECTOR 
C LEVEL. 

114 TOPX=X ( K) 

TOPY=Y ( K ) 

TOPW=W(K) 

TOPT =T ( K ) 


20 CONTINUE 

COMMENT IF (PTHI.NG) WAS NOT ASSIGNFD A VALUE FOR THIS PARTICLE THEN. 

C INCREMENT (EXITPT) AND CONTINUE. IF A VALUE WAS ASSIGNED, 

C FIND THE ACTUAL PHYSICAL PATH LENGTH, INCREMENT (TOTC), AND 

C TEST TO SEE IF THIS WAS THF I ONGEST PATH LENGTH SO FAR. 

115 I F (PTHI.NG. EO. 0. )GO TO 19 
PTHLNG=PTHLNG*H(KS)/CTHETA 
123 01 =LPPTH*PTHI.MG+1 

TOTC (POS, 01 )=TOTC(POS,OI )+l 
I F ( PTHLNG . GT . LAPTH ) LAPTH=PTHLNG 
IS EXITPT(KS,POS)=EXITPT(KS, POS ) + 1 
GO TO 117 

COMMENT IF INSUFFICIENT TIME REMAINED WHEN (REMTIM) WAS CALLED, SET 
C (NPAR) TO THE ACTUAL NUMBER OF PARTICLES PROCESSED AND SET 

C THE FLAG, (TIME), TO ONE. 

106 NPAR=J-1 
TIME=1 

WRITE(6, 1005 ) NPAR 

COMMENT WRITE OUT THF. LONGEST POSSIRLE PATH LENGTH, THF LONGEST ACTUAL 
C PATH LENGTH, AND THE AVERAGE PATH LENGTH. 

117 WLPPTH=1 00 . /LPPTH 

AVEPTH=TOT PTH/EX I TPT ( NPI ANE , 1 ) 

WR I TE ( 6, 1006)WLPPTH, LAPTH, AVEPTH 
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COMMENT IF THE DETECTOR HAS THE FORM OF A RECTANGULAR PARALLEP I PED, 

C COMPUTE THE GEOMETRICAL FACTOR ANALYTICALLY AND WRITE IT 

C OUT . 

DO 22 J-l, NPLANE 

IF((T(J) .NE.TO) .OR. (W(J) .NE.WO))GO TO 118 

22 H0=H( J )+H0 
SOX=SORT(HO*HO+WO*WO) 

$o"y«sqrt (H0*H0+T0*T0) 

CGF=2*TO*SOX*ATAN CTO/SOX ) +2*WO*$OY*ATAN (WO/SOY )-2*HO*TO* 

AATAN (TO/HO) -2 *W0*H0*ATAN( WO/HO )-HO*HO*ALOG( (H0*H0* (H0*H0+ 
BT0*T0+W0*W0) ) / (SOX*SOX*SQY*SQY) ) 

WRITE (6/ 1007)CGF 

COMMENT FIND THE GEOMETRICAL FACTOR AND WRITE IT OUT. 

118 GF-3.141593*T0*W0*EX I TPT (NPLANE, 1)/NPAR 
WRITE(6, 1008)GF 

COMMENT WRITE OUT THE NUMBER OF PARTICLES THAT LEFT THE DETECTOR, BY 
C PLANE, AT EACH LEVEL. 

DO 23 J B 1, NPLANE 

23 WRITE(6,1009)J,(EXITPT(J,K),K-1,3) 

WR I TE ( 6, 1010 ) 

COMMENT FIND THE EFFECTIVE GEOMETRICAL FACTOR FOR ALL PARTICLES THAT 
C HAVE A PATH LENGTH OF (PBIN) OR GREATER. 

DO 24 J-1,100 

24 TOTC(4,J)»TOTC(l,J)+TOTC(2,J)+TOTC(3,J) 

SUMC (100)-TOTC(4,100) 

DO 25 J-l, 99 

25 SUMC (100-J)=TOTC(4, 100-J)+SUMC ( 101-J ) 

DO 26 J-1,100 

EFFGF-SUMC ( J ) *3 . 141593*T0*W0/NPAR 
IFCEFFGF.EQ.O. )GO TO 119 
PBIN-J/LPPTH 

COMMENT FILL (GPB I ,GTOT,GEFF) , WRITE OUT THE PATH LENGTH BIN, (PBIN). 

C THE NUMBER OF PARTICLES HAVING (PBIN) PATH LENGTH ACCORDING TO 

C WHERE THEY EXIT THE LAST DETECTOR LEVEL AND THE EFFECTIVE 

C GEOMETRICAL FACTOR FOR PARTICLES OF (PBIN) PATH LENGTH OR 

C GREATER. 

GPB I (J)-PBIN 
GTOT (J)=TOTC(4, J) 

GEFF( J )=EFFGF 

26 WRITE(6,1011)PBIN,TOTC(2,J),TOTC(3,J) ,TOTC (l,J),TOTC(4,J),EFFGF 
COMMENT PLOT (EFFGF,TOTC(4, *) ) VS. (PBIN). 

119 CALL G( 5, SUBH1,GPB I ,GEFF, J ) 

CALL G( 5, SUBH2, GPB I ,GTOT, J ) 

WR I TE ( 6, 1012) 

ZERO=NPAR 

COMMENT FIND THE COUNTS PER SQUARE DEGREE, (CPSD), AT EACH ANGULAR 
C BIN, (POLAR). 

DO 27 J= 1, 30 
QS=QS+AND( J ) 

POLAR=3*J 
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COMMENT FILL (GAND,GPOL,GCPS) FOR PLOTTING. 

GAND ( J )=AND(J ) 

GPOL ( J ) =POLAR 

CPSD=AND(J)/( 20626. 48*COS( ( J-l)*5. 235988E-2)-COS(J*5. 235988E-2) ) 
GCPS ( J )=CPSD 

COMMENT WRITE OUT THE ANGULAR BIN, THE NUMBER OF PARTICLES IN THAT BIN 
C AND THE (CPSD) , SUPPRESS ZEROS BY TESTING WITH (QS), SAVE THE 

C VALUE OF (J). 

WRITE(G, 1013)POLAR,AND(J),CPSD 
JS=J 

IF(QS.EQ.EXITPT (N PLANE , 1 ) )GO TO 120 
27 CONTINUE 

COMMENT PLOT (AND, CPSD) VS. (POLAR) 

120 CALL G(5, SUBH3, GPOL, GAND, JS) 

CALL G(5, SUBH4, GPOL, GCPS, JS) 

COMMENT WRITE OUT (XDIS,YDIS) AT EACH PLANE WHERE THEY WERE ACCUMULATED 
C BY THF BIN FROM THEIR RESPECTIVE AXIS, ( XB, YB ) . SUPRESS ZEROS. 

DO 28 .1 = 1,11 

I F( ( ( J . FQ, 1) , AND . (GRAPH ( 2 ) . EG . 0 ) ) . OR . ( ( J . NE . 1 ) . AND . (GRAPH ( J ) . EO . 0) 
A) )GO TO 10 
I F( J , EQ. 1 )GO TO 121 
I 3=GRAPH(J) 

I 2=ORAPH( J-D + l 
DO 29 1 1=1 2, I 3 

29 ZERO=7FRO-F.XITPT( 1 1,2)-EXITPT( 11,3) 

121 V/R I TF f G , 1 014 ) GRAPH ( J ), ZERO 
DO 30 K= 1, 50 
XR=K*TGR(J)/100. 

YB=K*WGR(J)/inn. 

COMMENT FILL (GXB, GYB, GXP I , GYD I ) FOR PLOTTING. 

GXB ( K ) =XB 
GYB(K)=YB 
GXD I (K)=XDIS(K, J) 

GYD I ( K) =YD 1 S ( K, J ) 

30 WRITE(0, 1015)K,XB,XDIS(K, J),YB,YDIS(K, J) 

COMMENT PLOT (XDIS(*,.D) VS. (XB) AMD (YDIS(*,J)) VS. (YB) 

CALL G ( 5, SUBH5, GXB, GXD 1,50) 

CALL G( 5, SUBHG, GYR,GYD 1,50) 

WRITE(G, lOlG)GRAPH(J) 

TTANG=0 

COMMENT IN THIS LOOP FIND THE TOTAL NUMBER OF PART I CLFS, (TAMG), THAT 

C WERF IN EACH ANGULAR RIM, (ANGLE), AT EACH PLANE LISTED IN 

C (GRAPH ( K) ) AND FILL (GANG,GTAN) FOR PIOTTING. ZEROS ARE 

C SUPPRESSED WHEN WRITING OUT (ANGLE, ANGD I S, TANG) . . 

DO 31 K= 1, 30 
TANG=0 
ANGLE =3 . *K 
GANG(K)=AMGLE 
DO 32 1=1,4 
DO 32 M= 1, 4 
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32 TANG-TANG+ANGO I S (M, L, K, J ) 

WRITE (6, 1017) ANGLE, ( ( ANGD I S (M, L, K, J ) , L=l, 4) , M=l, 4 ) ,TAMG 

TTANG=TTANG+TANG 

GTAN(K)=TANG 

KS = K 

I F (TTANG. FO. ZERO)GO TO 122 
31 CONTINUE 

COMMENT FIND THF SIZE OF THE ARAY TO BF. TRANSM I TED INTO (G), THEN PLOT 
C (TANG) VS. (ANGLE). 

122 S I ZE=ANGLE/3 . 

CALL G ( 5, SUBH7, GANG, GTAN, SIZE) 

COMMENT ^^I^TH IS LOOP WRITE OUT (ANGDIS) AS (ANGD1, ANGD2 ) IN 9-DEGRF.E 
C RING INSTEAD OF THF USUAL 3-DFGREF BINS, IN A MATRIX FASHION, 


SUPPRESSING ZEROS. 

DO 33 K=l,10,2 
I F ( 3* ( K-l) ,GF. . KS )GO TO 28 
DFG1=K*9. 

DFG2“(K+1)*9. 

V/R I TE ( 6, 1019 )DFG1, DEG2 
DO 33 L=l, 4 
DO 34 M»l,4 
ANGDl(M) B 0 
ANGD2(M) =0 
DO 35 M-1,4 
DO 35 N-1,3 

ANGD1(M)=ANGD1(M)+ANGD IS(M,L,N+(K-1)*3,J) 

ANGD2 (M) =ANGD2 (M)+ANGD IS(M,L,N+(K-1)*3+3,J) 

WRITF(C, 1020)L, (ANGDl(M) ,M-1, 4) , L, (ANGD2 (M) ,M=1, 4 ) 

CONTINUE 
CONTINUE 
STOP 

FORMAT ( 1 OTHE TIME IN MIN. REQUIRED FOR THIS DATA SET WAS',F6.2, 
A' CPU ' , F4 . 2, ' I/O') 

1001 FORMAT ( 20A4/I 2,2F10. 4, 110, 15,1012) 

1002 FORMATC1DATA SET NO. ', I5//1X,20A4/'0NUMBER OF PLANES 


34 


35 

33 

28 

10 

102 

1000 


AST 


OF PARTICLES = IIO/'OIY = 
H( I ) 0 ' , 2( 2X, F10. 4) ) 


I 5/ * 0 I = 


FIR 
T ( I ) 


12/' ONUMBER 
B ® U(l) - 

1003 FORMAT (3F10. 4 ) 

1004 FORMAT ( I4,3(2X,F10.4)) 

1005 FORMAT ( ’OTHERE IS NOT ENOUGH TIME 
A NUMBER OF PARTICLES = ',110) 

1006 FORMAT( ' OLONGEST POSSIBLE PATH » ' , F 10 . 4/ ' OLONGEST ACTUAL PATH = 
A, F 10 . 4/ ' OAVE . PATH = ',F10.4) 

1007 FORMAT( 'OCOMPUTED GEOMETRICAL FACTOR = ',1PE12.C) 

1008 FORMAT( ' OGEOMETR I CAL FACTOR = ' , 1PE 12 . 6/ ' ODETECTOR 

ADE FRONT') 

1009 FORMAT (4X,I2,5X,3(I7,2X)) 

1010 FORMAT Cl PATH-LENGTH SI DECTS FRTCTS 
A'EFFGF') 


FOR MORE PARTICLES OR DATA 


BOTTOM 


SI 


TOTCTS ' , 4X, 
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1011 FORMAT ( 1H , F8 . 2, 4 I 10, F 13 . 2 ) 

1012 FORMAT ( ' 1POLAR ANGLE COUNTS OUT BOTTOM COUNTS PER SO-DEG') 

1013 FORMAT ( 111 ,I7,14X,16,12X,1PE12.6) 

1014 FORMAT ('101 STURP>UT I ON OF INCIDENT PARTICLFS AT PLANE NO. ',12/ 
A'OTOTAL PARTICLES AT THIS PLANE = ',110/' OB IN X BIN NO. OF PAR 
BTITCI.ES Y BIN NO. OF PARTICLES') 

1015 FORMAT (I4,2(2X,F8.2,4X, I 10, 3X ) > 

1016 FORMAT ( ' 1ANGULAR DISTURBUTION OF INCIDENT PARTICLES AT PLANE NO. ' 

A, 1 2/ ' 0 ',51X, 'CORDINATES (X, Y) '/' ANGLE (1,1) (1,2) (1,3) 

B( 1, 4 ) (2,1) (2,2) (2,3) (2,4) (3,1) (3,2) (3,3) (3,4) (4, 

Cl) (4,2) (4,3) (4,4) TOTAL') 

1017 FORMAT(2H , F4 . 1, 16 ( IX, I 6 ) , 2X, I 8 ) 

1018 FORMAT ( '1ANGULAR DISTURBUTION AT 9. DEG BINS') 

1019 FORMAT ( 1110, 2 ( 14X, F3 . 0 , ' DEG',32X)/1H ,2('Y/X 1 2 

A 3 4 ' , 16X ) /3H ,38('*'),15X,38('*')) 

1020 FORMAT ( 111 ,11,' * ' , 4 ( IX, I 7, IX) , ' * ' , 13X, 1 1, ' * ' , 4 ( IX, I 7, IX ) , ' * ' ) 

END 
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•SUBROUT I N F. G(TYPE,HFAD,A,B,5IZE) 

REAL* 8 HEAD ( 5 ) 

LOG I CAL*1 CWR( 103 ) 

REAL A RAY ( 101) , BARAY ( 101) , A(S I ZE ) , BCG I 7.F ) 

INTEGER ! ARAY (101), I GARAY ( 101 ) , Q.W, 0( 11 ) / TYPE , V/A, EXPB, AD I 


F,BDIF, 


ATYP, S I ZE, H I ST 
COMMENT IF (HIST) 
IF (HIST) 
GRAPHS. 

H I ST=1 


FOIIALS TO ONF THEN THIS SUBROUTINE MAKES HISTOGRAMS, 
DOSE NOT EQUAL TO ONE THEN THIS SUBROUT I NE ..MAKES 


.1 = 101 
AMAX=0. 


BMAX=0 . 

DO 10 1=1,101 
A RAY ( I )=0. 00001 
BARAY ( I ) =0 . 00001 

I F( ( I . LE. SIZE) .AND. (A( I ) . GT . . 00001 ) ) ARAY ( I ) =A ( I ) 

IF((I. LE. SIZE). AND. (B(I).GT.. 00001)) BARAY ( I ) =R ( I ) 

10 CONTINUE 

DO 11 1=1,101 

!F( (ARAY (101- I ) .LF.. .00001) .OR. ( BARAY ( 101- I ) . LE . . 00001) )GO TO 100 

11 J-101-1 

100 DO 12 1=1, J 

I F (ARAY ( I ) . GT . AMAX ) AMAX=ARAY ( I ) 

I F (BARAY ( I ) . GT . BMAX ) BMAX=RARAY ( I ) 

12 CONTINUE 

I F ( ( AMAX .IE.5.).OR.( RMAX . LF . 5 . ) ) GO TO 101 
WR I TE ( 6, 1000 )HEAD 
I F (TYPE . LT . 5 )GO TO 106 
TYP =1 

I F (BMAX . GE . 1 . E5 )TYP =TYP +2 
I F ( AMAX . GE . 1 . E5 )TYP =TYP +1 
106 GO TO (102, 103, 104, 105), TYP 

102 TYP=1 

AD I F=AMAX/100+ .995 
BD I F= RMAX/ 100+ .995 
CWR( 102 ) = 64 
CWR( 103) =64 
EXPP=1 

UR I TE ( 6, 1001 ) 

DO 13 1 = 1,. I 

I ARAY ( I ) = ( ARAY ( I )+ADIF/2. )/ADIF+l. 

13 I RARAY ( I )= (BARAY ( I )+RD I F/2 . )/BD I F+l 
GO TO 107 

103 TYP=2 
AD I F = 1 

RD I F=RMAX/100+ . 995 
CUR( 102) =64 
EXPB= 1 

CUR( 103)=197 
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WR I TE ( 6, 1002 ) 
no 14 1=1, J 

I ARAY ( I )= (ALOGIO (ARAY ( I ))+.05)*10+l 

14 I BARAY ( I )= (BARAY ( I ) + BD I F/2 . ) /BD I F+l 
GO TO 107 

104 TYP=3 

AD I F=AMAX/100+ .995 
BD I F= 1 

CURC102)=197 
CWR( 103 )=64 
EXPB=10 
WRITE(G, 1003) 

DO 15 1=1, J 

I ARAY ( I ) =(ARAY ( I )+AD I F/2 . )/AD I F+l. 

15 I BARAY ( I )=(AL0G10(BARAY( I ) )+.05)*10+l 
GO TO 107 

105 TYP=4 
AD I F = 1 
BD I F=1 

CWR(102)=197 
CWR( 103 )=197 
EXPB =10 
WR I TE ( 6, 1004 ) 

DO 10 1=1, J 

I ARAY ( I ) = (ALOGIO (ARAY ( I ) )+ . 05) *10+1 

16 I BARAY ( I )= (ALOGIO (BARAY ( I ) )+ . 05) *10+1 
107 I F (H I ST.NE . l)GO TO 115 

DO 22 1=1,101 
22 CWR( I )=64 
115 DO 17 1=1,101 

I F (H I ST. EQ. l)GO TO 119 
DO 18 d-1,101 

18 CWR(J)»64 

119 WA =BD I F* ( 101- I ) / F.XPB 
DO 19 K-1,11 

IF(I.EQ.10*(K-1>+1) GO TO 113 

19 CONTINUE 
WRiTE(6, 1005) 

GO TO 114 

113 WRITE (6, 1006 )CWR( 102 ) , WA 

114 DO 20 J-1,101 
QW=0 

I F ( I BARAY (J ) . EQ . 102- I )GW= I ARAY ( J ) 

I F (OW. GT . 0 )CWR( OW) =92 

20 CONTINUE 

17 WRITE (6, 1007) (CWR(J ),J=1,101) 

DO 21 1=1,11 

IF(TYP/2.EQ.TYP/2. )GO TO 116 
Q( I )=ADI F*( I -1)*10 
GO TO 21 
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lie Q( I ) = I -1 
21 CONTINUE 

WRITE (G,1008)(CWR(103),Q(I), 1=1, ID 
101 RETURN 

1000 FORMAT ( 1H1, 5A8) 

1001 FORMAT ('OTHIS IS A SCALER PLOT') 

1002 FORMAT ( ' OTH IS IS A SEMI -LOG PLOT IN X' 

1003 FORMAT ('OTHIS IS A SEMI -LOG PLOT IN Y 

1004 FORMAT ( ' OTU I S IS A LOG-LOG PLOT') 

1005 FORMAT(9X, ' I ' ,10<9X,'I')) 

1006 FORMAT ( 1H , Al, I 7, 101 ( ' - ' ) /10H+ 

1007 FORMAT ( 9H+ , 101A1) 

1008 FORMAT ( 1H , 5X, 11(A1, 1 3, 6X) ) 

END 


) 

') 

M0(9X, ' I ')) 
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DO 21 
L = 2, 11 


GRAPH (L) 
. EQ. K 


QX = ABS(XiK) ) • 100./T(K) + 1 
QY = ABS (Y (K) ) * 100. / W (K) + 1 
XDIS (QX, LS) = XDIS (QX, LS) + 1 
YDIS (QY. LS) = YDIS (QY, LS) + 1 
ANGX « ABS (X (K) * 8 / T (K) + 1 
ANGY = ABS (Y (K) * 8 / W (K) + 1 
ANGDIS (ANGX. ANGY. ANGT, LS) = 
ANGDIS (ANGX, ANGY. ANGT, LS) + 1 


PTHLNG = H (K) / CTHETA 
TOTPTH *= TOTPTH + PTHLNG 
QI = THETA ' 19.09859+ 1 
AND (QI) = AND (QI) + 1 


61 


64 












































XB - K * TGR |J) / 100. 
YB » K * WGR (J) / 100. 
GXB(K) - XB 


GYB(K) - YB 
GXDI (K) » XDIS(K.J) 
GYDI <K) - Y0IS (K, J) 

— 

K, XB.XDIS(K.J), 
YB. YDIS(K.J) 



TANG - 0 
ANGLE - 3. * K 
GANG (K) - ANGLE 



ANGLE, UANGDIS 
(M, L, K, J), L ■ 1, 4), 
M - 1. 4), TANG 



TTANG = TTANG + TANG 
GTAN (K) » TANG 
KS * K 
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STOP 


C-ll 



